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Abstract 

We consider statistics of the disruption and Lyapunov times in an hierarchical 
restricted three-body problem. We show that at the edge of disruption the 
orbital periods and the size of the orbit of the escaping body exhibit Levy 
flights. Due to them, the time decay of the survival probability is heavy- 
tailed with the power-law index equal to —2/3, while the relation between the 
Lyapunov and disruption times is quasilinear. Applicability of these results 
in an "hierarchical resonant scattering" setting for a three-body interaction 
is discussed. 
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1 Introduction 



Notwithstanding a three-centennial progress in the studies of the three-body 
problem, disruption of a three-body gravitational system still remains an 
enigmatic dynamical process. The statistics of the disruption times is of spe- 
cial interest. In particular, it is interesting whether the computed Lyapunov 
exponents of the motion are related in some way to the disruption times or, 
more generally, to the times of sudden changes in the orbital behavior. (The 
inverse of the maximum Lyapunov exponent, — the Lyapunov time Tl, gives 
the timescale of exponential divergence of close trajectories in phase space; 
it characterizes the level of predictability of the motion. The disruption time 
is the system lifetime as a bound system.) 

Much numerical-experimental work has been done on this subject in the 
last two decades. A prominent relationship in the statistics of sudden or- 
bital changes in gravitational systems has been found and confirmed. This 
relationship consists in the close-to-quadratic character of the dependence 
of the time T r of a sudden orbital change on the Lyapunov time Tl. (The 
time of a sudden orbital change is dynamically equivalent to the recurrence 
time for entering some domain in the phase space of a Hamiltonian system, 
therefore it is designated T r here and in the following.) On the basis of 
extensive numerical experiments in the dynamics of three-body and more 
complicated gravitational systems (in particular, in asteroidal dynamics), it 
was argued in [TJ [2j [3] that the time T r could be statistically predicted by 
means of computation of the maximum Lyapunov exponent. A simple "uni- 
versal" statistical relationship was established to exist between the computed 
time T r and the computed Lyapunov time T L (the inverse of the computed 
maximum Lyapunov exponent): 

T r oc Tl (1) 

with (3 = 1.7-1.8 typically, though considerable dispersion of the statistical 
data was usually present. The same kind of relationship with /3 w 1.9 was 
found in [I] in simulations of the dynamics of the Kuiper belt of asteroids. 
Similar dependences were found in [3 |6] for the chaotic behavior of aster- 
oids moving close to the mean motion resonances 2/1 and 3/1 with Jupiter. 
Tsiganis et al. [7j studied chaotic diffusion and effective stability of Jupiter 
Trojans, and in a large chaotic region, surrounding the stability zone in the 
motion of asteroids in the Trojan swarms, obtained a statistical relationship 



2 



between the escape and Lyapunov times, similar to the "universal" close-to- 
quadratic one. 

However, recent numerical-experimental statistical studies [319] of the dis- 
ruption process in three-body systems revealed a new — quasilinear — kind 
of relation between the Lyapunov times and the disruption times. Mikkola 
and Tanikawa [8J explored correlation of the Lyapunov times Tl and the 
disruption times T^ in the equal-mass three-body problem. In extensive nu- 
merical experiments they computed the disruption and Lyapunov times for 
the motion with randomized initial conditions, and showed that the system 
lifetimes, as a bound triple, and the Lyapunov times were correlated, the 
"Tl — Td" relationship at long lifetimes being close to linear. Urminsky and 
Heggie |Sj explored correlation between the Lyapunov and disruption times 
in an hierarchical three-body problem in a setting of the three-body problem 
different from that used in [8], but obtained similar results. Namely, they 
considered the Sitnikov problem [10] and, as an outcome of extensive numer- 
ical experiments, obtained a two-part power-law relationship with the second 
part (tail) close to linear one, similar to that found in [5] for the equal-mass 
three-body problem. 

Though the results [HI [9] on the "Tl — Td" relationship are similar, they 
diverge significantly from the earlier results JTJ [2j [3j HJ EE El [7] obtained for 
the statistics of sudden orbital changes. 

In this paper we identify causes that lead to the quasilinear "T L — T d " 
relationship in the three-body problem at the edge of disruption, as opposed 
to the earlier found close-to-quadratic "Tl — T r " dependence in resonant 
dynamics, and explain its nature as an effect of Hamiltonian intermittency. 



2 Two kinds of Hamiltonian intermittency 

The close-to-quadratic relationship ([1]) between the recurrence and Lyapunov 
times in general Hamiltonian systems with divided phase space (i.e., with 
regular and chaotic kinds of motion both present) has a straightforward ex- 
planation if one takes into account the so-called "stickiness phenomenon" 
[61 (TTJ: the correlation arises because the trajectories sticking to chaos bor- 
der have large recurrence times (the times of return to some domain of phase 
space), and, since they mimic regular ones, they have also large computed 
("local") Lyapunov times. The relationship between the recurrence time T r 
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and the local Lyapunov time Tl on each recurrence was shown to be close to 
quadratic [5J: 

T r oc T L 2 , (2) 

i.e., (3 m 2 in Eq. ([1]). This relationship was derived in [6J for the motion 
inside the chaotic layer near the separatrices of nonlinear resonance in the 
perturbed pendulum model, and this indeed covers a lot of possible applica- 
tions, including the chaotic motion of asteroids near mean motion resonances. 
The separatrix map theory [121 EH EH] was used for deriving relation (j2J). 

Generally speaking, the emergence of the algebraic "Tl — T r " relation- 
ship is due to Hamiltonian intermittency. The phenomenon of intermittency, 
i.e., chaotic behavior intermittently interchanged with close-to-regular one, 
is well known in studies of dissipative dynamical systems [2]. Much less it 
is discussed in connection with Hamiltonian dynamics. There are two kinds 
of Hamiltonian intermittency known. The first one [15] takes place when the 
motion is "stochastized" at encounters with the separatrix, while most of the 
time (far from the separatrix) the motion looks regular. The separatrix in 
the simple model [T3], as well as in the map models (JU E]) considered in this 
paper, is just the singular line y = separating bound and unbound types 
of motion. The first kind Hamiltonian intermittency is inherent to adiabatic 
chaos; see [161 E] and references therein. The second kind Hamiltonian inter- 
mittency [18] takes place when a fractal "chaos-order" boundary is present 
in phase space. Sticking to the border results in long time segments of close- 
to-regular behavior; they are interrupted by prominently chaotic dynamical 
events, when the trajectory leaves the border neighborhood [HI [201 El 118] . 
Thus the physical cause for the Hamiltonian intermittency of the first kind 
consists in close encounters of a trajectory with slowly pulsating separatrix, 
while the physical cause for that of the second kind consists in sticking of a 
trajectory to chaos border. Statistical properties of these two kinds of phe- 
nomena are very different, as demonstrated further in this paper. This is 
just the second kind Hamiltonian intermittency that is responsible for the 
quadratic "Tl — T r " relationship fl2]) [6J. In what follows we show that in 
certain circumstances a quasilinear "Tl — T r " relationship can arise, and 
this is due to the first kind Hamiltonian intermittency. 
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3 General separatrix maps 



The nonlinear pendulum provides a model of nonlinear resonance under def- 
inite conditions [T2| IT4"] . The motion in the vicinity of the separatrix of 
the nonlinear pendulum (nonlinear resonance) is described by the separatrix 
map (121 [H] . We write it in the form adopted in [6] : 

Vi+i = Vi + sinxi, 

x i+ i = Xi + A In \y i+1 \ + c, (3) 

where y denotes the normalized relative pendulum's energy, x is normalized 
time, the constants A and c are parameters. 

Consider a map similar to the separatrix map ([3]), but with a power- law 
phase increment instead of the logarithmic one: 

y i+1 = y { + sin Xi, 

= Xi + A|y i+ i|" 7 , (4) 

or, in an equivalent form usually used, 

w i+ i = Wi + W sin Ti, 

Ti+l = Ti + K\lU i+1 \^. (5) 

Map (jSJ) has two parameters, W and k, instead of one parameter A in map (j4]); 
apart from the 7 parameter. The two-parameter map (jSJ) is reduced to the 
one-parameter map (jlj) with A = kW~^ by means of straightforward substi- 
tution w = Wy, t = x. 

A number of mechanical and physical models are described by maps (jlj) 
and ((Sj) with rational values of 7. The values of 7 = 1/4 and 1/3 correspond 
to the Markeev maps [211 122] for the motion near separatrices of resonances 
in two degenerate cases; 7 = 1/2 gives the "L-map" [12] for the motion of 
a non-relativistic particle in the field of a wave packet, this value of 7 also 
gives a map for the classical Morse oscillator driven by time-periodic force 
[23] ; 7 = 1 gives the Fermi map [2HQ3] for the Fermi acceleration mechanism 
for cosmic rays; 7 = 3/2 gives the Kepler map [25l [26J [271 (28J [23] f° r a number 
of physical and astronomical applications; 7 = 2 gives the "ultrarelativistic 
map" P2] for the motion of a relativistic particle in the field of a wave packet. 
All these maps describe dynamical behavior in the vicinities of separatrices 
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of corresponding models. In the case of the Kepler map the separatrix (the 
line y = 0) separates the bound and unbound states of motion. 

The Kepler map (Eqs. (jl]) or (j5]) with 7 = 3/2) was derived and analyzed 
in [25l 126] 1271 128] in order to describe the chaotic motion of the Halley 
comet and, generally, the motion of comets in highly eccentric orbits. The 
motion model for the Kepler map consists in the assumption that the main 
perturbing effect of planets (Jupiter first of all) is concentrated when the 
comet is close to the perihelion of its orbit. The y variable has the meaning 
of the normalized orbital energy of the comet, while x is the normalized time. 
One iteration of the map corresponds to one orbital revolution of the comet. 
This means that the map time unit, corresponding to one iteration, is not 
constant. The increment of real time per iteration is Axj = Xi — Xi-\. 

The Kepler map is known to describe dynamics in several different set- 
tings of an hierarchical restricted three-body problem: in the external re- 
stricted planar [25] 126] 128] and strongly non-planar [29] problems in cometary 
dynamics; also in the abstract Sitnikov problem, where the tertiary moves 
along the perpendicular to the orbital plane of the main binary. The validity 
of the map in the last case follows from the work by Urminsky and Heggie 
[9], who considered a variant of the Sitnikov problem and derived a map 
(see Eqs. (11) in [9]), describing the dynamics in this problem. This map 
is straightforwardly reducible to the Kepler map (jl]); in particular, using 
formulas given in [9] for the parameters of their map, it is easy to derive a 
formula for the parameter A in Eqs. 01]): 

A = 27r(2.029e)- 3/2 « 2.17e~ 3/2 , (6) 

where e is the eccentricity of the central binary. So, for the eccentricity e 
equal to 0.1 and 0.6 (chosen as representative in [9]) one has A ~ 70 and 
5, respectively. We choose these two values of A for illustrative numerical 
experiments in what follows. 

4 Levy flights at the edge of escape: 
the distribution 

The Levy flights, i.e., the increments (in a kind of random walk) that have a 
heavy-tailed distribution, is a well-studied subject with various applications. 
In Hamiltonian dynamics, they were thoroughly considered in connection 
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with what we call the second kind Hamiltonian intermittency; see, e.g, [3U] 
and references therein. In celestial mechanics, Levy random walks arising due 
to close encounters of bodies were considered in the orbital energy evolution 
of comets [3T] . 

Usually, Levy flights are considered in random walks with steps possible 
in at least two (forward and back) directions; but the Levy distributions 
considered below are one-sided: the increments (the orbital periods of the 
escaping body and the recurrence times) are positive. In a general statistical 
setting, one-sided Levy flights were considered in [32], [33] , where exact results 
for the first passage time and leapover statistics were obtained. 

Both kinds of Levy flights explored below in the framework of the three- 
body problem are not due to encounters of bodies. One of these two kinds is 
due to encounters of a trajectory with the separatrix in phase space; hereafter 
we call such flights the "Levy flights of the first kind", or LF1, since this 
phenomenon corresponds to the Hamiltonian intermittency of the first kind. 
Levy flights of another kind, which we call the "Levy flights of the second 
kind", or LF2, are due to sticking of a trajectory to chaos border. In other 
words, these flights arise due to the second kind Hamiltonian intermittency. 
They effect the duration of long Poincare recurrences, but not (practically) 
the orbital periods. 

In the orbital behavior of the tertiary in the planar restricted three-body 
problem which we consider in what follows, the LF1 appear as sudden jumps 
in the orbital size and period, while LF2 appear as long sequences of orbital 
revolutions with almost constant orbital size and eccentricity. 

LF1 still take place when encounters of bodies are impossible. E.g., in 
the cometary dynamics model, described by the Kepler map, [25, 26J the 
perihelion distance of the comet can be greater than the semi-major axis of 
the orbit of the secondary by any amount, but encounters with the separatrix 
with or without crossing it can still take place (if the cometary orbit is 
chaotic) and, as we shall see below, shape the statistics of recurrences. 

In the dynamics of map (jl]), the LF1 and LF2 coexist (if A is large enough, 
this condition is considered below). In Fig. [TJ a fragment of a chaotic trajec- 
tory of map (jl]) with 7 = 3/2 is shown, demonstrating LF1 (narrow peaks 
in the variation of the orbital period P or b of the tertiary, on the left of the 
graph) and an LF2 (the oscillatory low "plateau" in the variation of y and 
P or b, on the right). When y hits close to zero, a jump in the orbital period is 
observed. When the trajectory sticks to chaos border, both the energy and 
the orbital period oscillate near some low constant value. 
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The stickiness effect determines the character of the distribution of Poincare 
recurrences on large timescales: it is algebraic [SUES]. The algebraic decay 
in the recurrence statistics in Hamiltonian systems with divided phase space 
was considered, in particular, in [3U ESJ [131 ESJ EI] , starting with a pioneering 
work by Chirikov and Shepelyansky [31]. Chirikov [13] . using his resonant 
theory of critical phenomena in Hamiltonian dynamics, justified a value of 
3/2 for the critical exponent a in the integral distribution 

F oc T~ a (7) 

of recurrences. The integral distribution F(T T ) is defined here as the relative 
share of the recurrences with the duration greater than T r in the whole sam- 
ple. In a recent paper [37] the algebraic decay of Poincare recurrences was 
explored statistically on the basis of large computational data on behavior of 
various Hamiltonian systems. These numerical experiments showed system- 
dependent power-law exponents, but the mean "universal" exponent turned 
out to be well-defined and equal to 1.57 ± 0.03, somewhat differing from the 
standard 3/2 value. Venegeroles [M] reports a value, equal to 1.54 ±0.07, re- 
sulting from averaging independent results of numerical studies of a number 
of various Hamiltonian systems; see Table 1 and references in [38]. 

In celestial mechanics, the algebraic decay was observed in numerical ex- 
periments on asteroidal dynamics [T$ll2"0"] . These experiments were performed 
in the framework of the restricted three-body problem Sun- Jupiter-asteroid. 
It was shown that the tail of the distribution of duration T r of intervals be- 
tween jumps of the orbital eccentricity of asteroids in the 3/1 mean motion 
resonance with Jupiter is algebraic: 

F oc T~ a (8) 

with a ~ 1.5-1.7 typically. This was interpreted in [20] as an effect of sticking 
of the chaotic orbits to chaos border in the divided phase space. In other 
words, this is an effect of the second kind Hamiltonian intermittency. 

On the other hand, Dones et al. [39J studied the escape times in the highly- 
eccentric chaotic cometary dynamics in the Solar system (the perturbations 
due to the four giant planets were taken into account) and reported on the 
algebraic tails of the integral distributions with the power-law index equal to 
0.8 ± 0.2. This behavior has not been theoretically interpreted up to now; 
we shall see that it can be straightforwardly interpreted as an effect of the 
first kind Hamiltonian intermittency. 
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The decay with a ~ 1.5 is expected when LF2 dominate over LF1 in 
the long recurrences. A different kind of the decay law was observed in the 
behavior of the Kepler map by Borgonovi et al. [40J. They analyzed long- 
time decay properties of a Kepler map describing a one-dimensional model 
of hydrogen atom in a microwave field, and by means of rigorous deduction 
found a T d _2//3 law for the first approximation for the time decay of the survival 
probability in the case of the escape times measured in real (constant) time 
units. This law was confirmed by them in computed statistics; see Fig. 1 
in [ID] . On the opposite, when the escape times were measured in map 
(fictitious) time units, the usual T d 3 ^ 2 law was observed; see Fig. 2 in |4"D] . 

Apart from the rigorous treatment [3D], several heuristic deductions of 
the T d 2//3 law are available in relevant problems [311 1321 133] . Hut [32] derived 
an heuristic T d 2 ^ 3 law as a lower bound for the time decay of the survival 
probability in a general "hierarchical resonant scattering" [331 132] setting 
for a three-body interaction (where the masses of "stars" are arbitrary). 
Malyshkin and Tremaine [33] derived the T d _2,/3 law for the time decay of the 
survival probability for cometary ensembles in the Solar system. Schlagheck 
and Buchleitner [31] derived the T d 2//3 law for the time decay of the survival 
probability in an autoionizing configuration of chaotic helium. 

In all these approaches, the T^ 2 law was derived without using any Ke- 
pler map. Two basic assumptions were always made explicitly or implicitly, 
(i) that the distribution of ejection energies is flat or smooth in the neigh- 
borhood of the energy threshold E = 0, (ii) that the asymptotic decay of the 
survival probability is the same as the tail of the distribution of the orbital 
periods of the escaping body. While the first assumption is reasonable (in 
view of the complete ergodicity of the motion near the threshold, see anal- 
ysis given below for the Kepler map), the second one is solely hypothetical. 
E.g., the argumentation presented in [33] is that the escaping comets "remain 
bound until their second perihelion passage, after which they will normally 
be ejected within a relatively short time". However, given that all results 
[311 1321 133] coincide with the result of the rigorous treatment in [40J , we can 
infer that the second assumption seems to be also valid. 

Let us, using the same two assumptions, find a law for the asymptotic 
distribution of the phase increments per iteration (which are the orbital pe- 
riods of the escaping body in case of 7 = 3/2) in the general case of arbitrary 
7 in map (j3D- The phase increment is 

P = Axi = Xi - = Xlyil" 1 (9) 
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in real time units. In case of 7 = 3/2 the real time unit is equal to the orbital 
period of the central binary, divided by 2n. 

So, \yi\ = \ l hp- l /i _ The motion in the close vicinity of the separatrix 
is locally ergodic, i.e., regular islands are absent. This follows from the fact 
that in the local (in the energy y) approximation of map (jlj) by the standard 
map the stochasticity parameter of the standard map tends to infinity when 
y tends to zero. The ergodicity implies that, for the y variable close to zero, 
the distribution function of y is flat: f(y) = const. One has 

const d|y| oc P'y^dP, (10) 
and the differential distribution function of P is 

/(P)ocP^- 1 . (11) 

The integral distribution is F(P) oc P ~ . 

For the Kepler map, 7 = 3/2 and P = P orD (the orbital period); so, the 
differential distribution is /(Port.) oc P orb and the integral distribution is 
P(-forb) oc -Porb'' 3 - From Kepler 's third law one has f(a) oc a~ 2 for the differ- 
ential distribution of the semi-major axis of the orbit; i.e., the distribution 
of the orbital size is also heavy-tailed, and the Levy flights are demonstrated 
in the process of disruption both in tertiary's orbital period and size. 

By assumption (ii), advocated above, law ( iTTj) coincides with the asymp- 
totics of the time decay of the survival probability. Alternatively, the same 
law follows from repeating the analytical treatment [10] for the case of arbi- 
trary 7 instead of the exponent 3/2 in the Kepler map: if one repeats the 

rigorous deduction [JDj (performed in jlU] for the Kepler map, 7 = 3/2) in 

_i 

the general case of arbitrary 7, one finds the distribution F{T r ) oc T r 7 for 
the recurrence times. Note that the basic assumption in this deduction is 
that the phases are randomized after each kick. 

Borgonovi et al. [10] explained the LF1 dominance over LF2 (in our terms) 
in real time statistics in the dynamics of the Kepler map, using argumentation 
based on the infinite measure of extended phase space near the separatrix. 
If one uses this argumentation in the general case of map (jlj) with arbitrary 
value of 7 > 0, it follows that LF1 should dominate for all 7 > 1 at least. 
However, while LF1 dominate indeed at 7 > 1, the transition to this domi- 
nation occurs at a value less than 1. Let us estimate this critical value. At 
7 < 7crit ; when LF2 dominate, the slope index is critical (in a different sense, 
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related to the "critical" structure at chaos border, see discussion above): 
Oi = a cr it ~ 3/2. At 7 > 7 CI . it , when LF1 dominate, the slope index a = I/7. 
These two curves 0(7) intersect at 7 = 7 cri t = l/a C rit ~ 2/3. 

What is the physical reason for the switch to take place at this point? 
For the recurrences forming LF2 with duration greater than T, the "total 
sojourn time" ~ TF(T) [35], i.e., ~ T~ aciit+1 . Analogously, in the LF1 case, 
the total sojourn time is ~ TF(T) oc T~~ +1 . LF1 asymptotically dominate, 

if the second sojourn time is greater than the first one: 

Hence the condition for the LF1 domination is 7 > l/a cr ; t . This is exactly 

what we have just derived for the point of intersection of the two curves 0(7). 

So, a critical non-trivial value 7 crit of the 7 parameter exists (if A is large 
enough, see below), such as the maps with 7 > 7 crit have LF1 dominating 
over LF2 in real time statistics, whereas at 7 < 7 cr i t LF2 dominate in both 
the real time and map time statistics. 

The quantity 7 cr i t equals 2/3 in the case of the standard value a = 3/2; 
and 7 cri t ~ 0.637 in the case of a = 1.57 computed in [37]. In particular, LF1 
dominate over LF2 in real time statistics in the long Poincare recurrences in 
the dynamics of the Fermi and Kepler maps, while for the ordinary separa- 
trix maps, Markeev maps and L-maps, if A is large enough, the tails of the 
recurrence distributions are LF2-dominated and their slopes do not depend 
on the choice of units (map units or real time units) in which the lengths of 
recurrences are measured. 

To illustrate the difference between distributions of recurrences measured 
in map and real time units, we present here examples of distributions ob- 
tained by iterating map (JIJ) with 7 = 3/2 (the Kepler map). In Fig. [2} the 
computed integral distributions F(T V ) of the recurrence times T r measured in 
map time units (iterations), are shown for A = 5 (bold line) and A = 70 (thin 
line). The number of iterations n- lt = 10 11 in both cases. The recurrences are 
counted at the line y = 0. The quantity F(T T ) is the fraction of recurrences 
longer than T r . The tail of the distribution in Fig. [2] follows approximately 
the power law with the slope index a ~ 1.5, as expected. It is similar, e.g., 
to the tails of distributions presented in Fig. 4 in [20] for the intervals be- 
tween eccentricity bursts of chaotic asteroidal trajectories in the 3/1 Jovian 
resonance, where the same critical dynamical mechanism (sticking to chaos 
border) is present. In Fig. [3J the same distributions as in Fig. [2] are shown, 
but the recurrence times are measured in real time units. The slope index 
for the tails is evidently equal to 2/3, as expected in this case. So, when real 
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time units are used, the tails of the distributions demonstrate a behavior 
that is different in slope and regularity, compared to the case of using the 
map time units. The reason is that LF1 dominate in the distribution tail in 
the first case, while LF2 dominate in the second case. 

As we have already mentioned above, the LF1 and LF2 coexist in the dy- 
namics of map (jl]), if A is large enough. This condition provides the existence 
of the global fractal chaos border and, consequently, the prominent sticking 
phenomenon. If A is small enough, map (jl]) can be reduced to a differen- 
tial equation describing regular trajectories, at all y far from the separatrix 
y — 0. This can be done analogously to the case of the ordinary separatrix 
map ([3]); the procedure is described in [TB]. Thus the motion far from the 
separatrix is locally regular, no global fractal chaos border exists; only LF1 
are possible. This is the realm of "adiabatic chaos" . 

What is the boundary value of A, separating the cases with and without 
global fractal chaos border? It can be estimated from the form of the A 
dependence of the maximum Lyapunov exponent L of map fl4]) . The L value 
increases with A, while A is small, but then saturates at some constant value 
(see [13 fig. 3] for the Kepler map case 7 = 3/2). The saturation takes place 
when the role of the global fractal chaos border becomes important in the 
dynamics. As follows from jl5l fig. 3], the boundary value of A for this map 
is ~ 2-3. In the case of the ordinary separatrix map, the transition value of 
A is ~ 0.5-1 (see (Ml fig. 3], also [TBI fig. 1]). By means of constructing the 
dependences L(X) for arbitrary < 7 < 2, it can be shown that the boundary 
value of A does not change much with 7 and, by the order of magnitude, is 
~ 1. Thus, if A ^> 1, we expect coexistence of LF1 and LF2 in the dynamics 
of map @, while if A 1, only LF1 are possible. 

Taking this into account (i.e., setting A to be large enough), one can ex- 
ploit the phenomenon of coexistence of LF1 and LF2 to estimate the value of 
the critical exponent a. Namely, one can compute the distribution functions 
(in real time) of the Poincare recurrences for map (J4]) for a set of values of 7 
with some step, and, on increasing 7, fix its transition value 7 cri t, when the 
tail of the integral distribution function starts to take the form characteristic 
for the LF1 statistics, i.e., the algebraic decay with the power- law index a 
equal to I/7. Then a = 1 / 7crit • Thus the procedure gives a new method for 
estimating the critical exponent a in Eq. ([7]). An example is given by the 
graphs in Fig. HI One can see that the realization of our method with only 
two values of 7, namely, 7 = 1/2 (corresponding to the L-map) and 7 = 1 
(corresponding to the Fermi map), clearly shows that 1 < a < 2. 
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Now let us investigate this effect in more detail by building the depen- 
dence 0(7) directly, i.e., by finding a numerically on a grid of values of 7. 
When LF2 dominate, this dependence is expected to be subject to deforma- 
tions arising due to marginal resonances at chaos border: at some intervals of 
7, the border of the chaotic layer and, consequently, the recurrence statistics 
are perturbed due to emergence of the marginal resonances; on the marginal 
resonances, see [IS] . In order to reduce maximally these perturbations, let 
us introduce a constant shift c in x in map (Jl]), in the following way: 

Vi+i = Vi + smxi, 

Xi+i = Xi + Myi+il^ + c. (12) 

The shift c is a new parameter, analogous to the parameter c in map At 
each value of 7 we adjust the value of c in such a way that the winding number 
of the last undestroyed global KAM curve at the border of the chaotic layer 
is approximately equal to the "golden mean" (5 1 / 2 — l)/2; in other words, it 
is maximally far from the main resonances. 

By means of linearization of map (jl]) or (fT2l in y it is straightforward to 
see that the value of y corresponding to the critical value of the stochasticity 
parameter K = Kq m 0.971635406 [12J[I11IIE] of the approximating standard 
map is ?/b = (7A/ 'Kq)^ 1 . This value corresponds to the border of the chaotic 
layer. The regularizing constant shift c is found by averaging the phase 
increment at the border of the chaotic layer; the resulting formula is 

c reg = n(5^ 2 - 1) - A( 7 A/^ G )-^ +1 ). (13) 

Since the form of the relationship at(y) is expected to be independent 
from A, when A is large enough, one can take any large value A ^> 1; we 
choose A = 10. We build the graph on the interval < 7 < 2 with the 
resolution A7 = 0.01, i.e., we make 200 measurements of a. At each value 
of 7 on the grid, the integral distribution built in logarithmic coordinates is 
linearly fitted on the interval 3 < log 10 T r < 6; T r is measured in real time 
units. The left border of the interval is chosen to be much greater than the 
point of transition from the initial behavior (which can be exponential or 
inverse square-root, see [36]) to the asymptotic one. A direct inspection of 
the constructed distributions shows that this transition occurs at log 10 T r < 2 
for all points on the 7 grid. The right border of the T r interval is defined to 
be much less than the values at which the distributions exhibit a drop due 
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to poor statistics at the very tails. At each point on the 7 grid, the number 
of iterations of the map is n it = 10 10 . 

In Fig. |5l the computed dependence 0(7) is shown for two cases, c = 
and c = c reg . The bold line corresponds to the map with c = c reg , whereas 
the thin line corresponds to the map with c = 0. One can see that in the 
first case the irregular perturbations of the dependence are indeed suppressed 
seriously, though they are not at all completely eliminated. In the both cases, 
the transition to the LF1 domination at 7 ~ 0.7 is evident, in agreement with 
our prediction. However, it is also evident that the numerical construction of 
the 0(7) dependence does not provide a high-precision tool for determining 
the o; cr i t value, because the perturbations of the theoretically "flat" behavior 
at 7 < 7 cri t are large, and, due to them, the transition point cannot be 
located precisely. What is more, the suppression of the border perturbations 
by adjustment of the c parameter might introduce systematic errors (yet 
unexplored theoretically) in the estimated value of a C rit- 

Another way seems to be much more promising here. To locate the tran- 
sition point, one can exploit the property of fluctuations of the recurrence 
time distribution in the case of the second kind Hamiltonian intermittency, 
as opposed to the perfectly unperturbed behavior in the case of its first kind 
(see Figs. HJ). In our fits of the T r distributions, the a values are determined 
with their standard errors, i.e., at each point of 7 one obtains a ± o as an 
estimate for the power-law index. In Fig. El high-resolution 7 dependences 
of the standard deviation a for the numerically determined power-law index 
a are shown for A = 5 (lower line) and A = 10 (upper line); c = c reg . The 
interval of 7 is taken to be in the neighborhood of the expected switch value; 
namely, 0.55 < 7 < 0.80. The resolution of the plot in 7 is A7 = 0.0005. 
All other details of the numerical procedure are the same as adopted for 
construction of Fig. [5J 

In the constructed dependences in Fig. [6j a transition from a linear decline 
from large values of a to a low-level horizontal plateau takes place in a narrow 
interval ~ 0.67 < 7 <~ 0.68 (corresponding to a cr i t ~ 1.5), in agreement 
with our prediction for the critical value of 7. Concluding, construction and 
analysis of the dependence 17(7) does indeed provide a perspective numerical 
tool for determining the a cr i t value. 
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5 Levy flights at the edge of escape: 
the "Tl — T r " relation 

As discussed in the Introduction, the power-law character of the dependence 
of the times of sudden orbital changes on the Lyapunov times [U El El lU 
El El E] was explained in [SI [EE] as a phenomenon of critical dynamics, i.e., 
an effect of motion near chaos border. This is an effect of the second kind 
Hamiltonian intermittency. 

Here we show that the first kind Hamiltonian intermittency in the three- 
body problem at the edge of disruption leads to a quasilinear "T L — T r " 
relation, when the times are measured in real time units. As follows from 
Eqs. @, the length of a Poincare recurrence in real time units is 

n n 

TM = £A^ = A£ N - 3/2 , (i4) 

i=l i=l 

where n = T r is the duration of the Poincare recurrence in map time units 
(iterations); the time increment Axi = Xi — and the iteration i = 1 is 
set for the start of the recurrence. 

Let us consider the finite-time maximum Lyapunov exponent (of an origi- 
nal dynamical system) calculated for a recurrence in the case when real time 
units are used. To obtain it, one should divide the calculated finite-time 
maximum Lyapunov exponent L of the map by the average (on the recur- 
rence) length (Pmap) = Tj m ' /T x of the map iterations in real time (see, e.g., 
P3|). Hereafter we denote (P^l) by q: 

T->(ru) 1 n \ n 

The maximum Lyapunov exponent referred to real time units is 

= -. (16) 
Q 

Equivalently, for the Lyapunov time one has 

Tt ] = qTl- (17) 

On the other hand, 

T r ( m ) = qT r . (18) 



15 



The quantity q = Tj™' /T r is the ratio of two random variables, distributed as 
power laws, as shown above. If the y value hits close to the separatrix, then, 
according to (Tl4l) . there is a jump in Tj TU \ but there is no jump in T r . If the 
magnitude of such jumps is much greater than the total range of the original 
"Tl — T r " graph in map time units (this is the case when LF1 dominate over 
LF2), the graphical relationship in real time units will be spread, due to the 
jumps, in the direction Tl = T r . So, when Tl are T r are expressed in real 
time units and LF1 dominate over LF2, the "Tl — T r " relationship becomes 
quasilinear. 

This inference is valid when LF1 dominate. For the ordinary separatrix 
map ([3]) and general maps (TjJ with 7 < 7 crit , where singularity is weaker, LF2 
dominate if A is large enough, and the generic relationship is not spread in the 
Tl = T r direction. This is the reason why the close-to-quadratic relationship 
[D 121 El lU El El CTJ , and not the quasilinear one, is present in the dynamics 
of minor Solar system bodies near resonances, where the ordinary separatrix 
map, but not the Kepler map, is relevant to the dynamics of interacting 
nonlinear resonances. 

To illustrate our theoretical inferences, let us compute the "T L — T r " 
relationships for the Kepler map in map time units (iterations) and, for com- 
parison, in real time units. The finite-time maximum Lyapunov exponent is 
computed for a recurrence. In Fig. [7J the computed relationship in the map 
time units is shown for A = 5, n- lt = 10 7 . One can see, that, judging by the 
general slope of the dependence in log-log scale, the dependence is far from 
being linear here. Its slope in log-log scale is much steeper: the power law 
index is equal to 1.5-2, as expected. Now let us measure the recurrence times 
in real time units. Apart from this change in the time units, the same (as in 
Fig. [7J dependence is built in Fig. El Now the relationship is evidently quasi- 
linear, in accord with our theoretical finding for the case of real time units. 
One can even see how the "spreading" mechanism operates: the dependence 
has a " V" form, where the left wing is much shorter and represents a remnant 
(left after spreading by LF1) of the generic close-to-quadratic relationship. 
Thus the general "composite" appearance of the "Tl — T r " diagram (in real 
time units) for the Kepler map mimics general structure of the "Tl — T r " di- 
agrams revealed in computations of the disruption process in the three-body 
problem in various settings; compare Fig. El with figures 2 and 3 in [8] or with 
figures 3 and 7 in [9] . 

The over-all appearance of the diagram in Fig. [7J looks rather irregular, in 
comparison with that in Fig. El The nature of this irregularity can be clarified 
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by means of construction of a spectrum of winding numbers [17]. It can be 
built here for the Kepler map analogously to construction of the spectrum 
for the ordinary separatrix map in [U] . The spectrum of winding numbers 
graphically demonstrates which of the resonant chains of islands produce the 
longest events of sticking. The winding number Q is formally defined for a 
recurrence as Q = Ax/n, where Aa; is the length of a Poincare recurrence 
(the duration of an interval between crossings of the separatrix), measured as 
the total sum of variations in the phase x taken modulo 2ir, and n = T r is the 
length of the Poincare recurrence measured in the map iterations. We build 
the spectrum of winding numbers by plotting log 10 n versus Q. In Fig. HI 
the spectrum of winding numbers computed for A = 5 and 77. it = 5 - 10 11 is 
shown. The Farey tree [IS] of resonances is evident. Consider some lowest 
order resonances m/n and m'/n' that are "neighboring", i.e., mn' — m'n = 1, 
then the lower level of the tree is made of "mediants" given by the formula 
m"/n" = (m + m')/(n + n'). In Fig. [91 the evident lowest order resonances 
are 1/7, 1/6 and 1/5. The mediants for them are 2/13 and 2/11. For the 
resonances 1/7 and 2/13 the mediant is 3/20; for the resonances 2/13 and 
1/6 the mediant is 3/19; and so on. All mentioned resonances produce visible 
peaks and are easily identified in Fig. |9j We see that the irregular structure 
in Fig. [7J is explained by overlapping of individual relationships for several 
sticky island chains. The overlap of these relationships at the given timescale 
of recurrences produces the observed irregularity. 



6 Range of applicability 

Let us consider in more detail the range of applicability of the presented 
results on the statistics of the disruption process. Clearly, these results are 
valid wherever the Kepler map 01]) description of the motion of the tertiary 
is valid. First of all, we assume that the orbit of the escaping body is highly 
eccentric, and its pericenter distance q is much greater than the size of the 
orbit of the main binary. 

Hence the inferred statistics of the disruption and Lyapunov times are ex- 
pected in an hierarchical restricted three-body problem, where the pericenter 
distance of the tertiary is much greater than the size of the orbit of the main 
binary. The eccentricity of the orbit of the main binary as well as the non- 
coplanarity of the three-body system do not play role due to the following 
reasons. As it has been already mentioned above, the Kepler map is known to 
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describe the highly eccentric motion of a tertiary in several different settings 
of an hierarchical restricted three-body problem: in the external restricted 
planar and strongly non-planar problems, also in the Sitnikov problem. If 
the pericenter distance is large enough, only one harmonic in the Fourier 
expansion of the energy increment is important and is enough to be taken 
into account [2"%| 139"]. 

A role of the mass parameter \x in the main binary has not been explored 
yet, but one may expect the validity of the Kepler map approximation for 
the motion again, when q is large, because the higher order harmonics in the 
energy increment expansion are exponentially small with q [28, 49J. More- 
over, the formula for the energy increment in Eqs. §5§ in the case of the 
circular-orbit main binary with the mass parameter \x = 1/2 (equal-mass 
binary) is similar to that in the case of the circular-orbit main binary with 
small /i. To clarify this point, let us compare the formulas. In the circular 
planar restricted three-body problem with small \i one has for the energy 
increment, if q ^> 1 (q is measured in the units of the semi-major axis of the 
main binary): 



in the case of retrograde orbits of the tertiary [281 formulas (3.16)]. In the 
case of the equal-mass main binary and the non-planar problem, it can be 
found from [50, formula (26)] in the restricted problem limit that 



where / is the inclination of tertiary's orbit. We see that the structure of the 
formulas is similar, putting aside some differences in the numerical values 
of the coefficients and the power-law indices. Note that formulas (fT9l |20~| 
l2~Tj) can be used, if necessary, to find numerical values of the A parameter in 
Eqs. (jlj), because these formulas provide estimates of the W parameter in 
Eqs. (jSJ), while the formula for k in Eqs. (jSJ) is trivial (it is simply the time 
normalization, see, e.g., [28]). 




(19) 



in the case of prograde orbits of the tertiary, and 




(20) 




(21) 
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The inferred statistics appear to be still valid in a more general "hierar- 
chical resonant scattering" [HJ E2] setting for a three-body interaction, where 
the masses of "stars" are arbitrary. Hut [12] derived an heuristic T^ 2 ^ 3 law 
as a lower bound for the time decay of the survival probability in this setting, 
and showed this law to describe well the tails of numerical-experimental dis- 
tributions. Note that this heavy-tailed distribution is in accord with an early 
finding by Agekian et al. [51] that the mean life-time of a general isolated 
three-body system is infinite. 

As we have seen above, two basic assumptions are necessary for heuristic 
derivation of the T d 2//3 law, namely, (i) that the distribution of ejection en- 
ergies is smooth in the neighborhood of the energy threshold E — 0, (ii) that 
the asymptotic decay of the survival probability is the same as the tail of the 
distribution of the orbital periods of the escaping body. These two assump- 
tions look rather plausible even in the general three-body problem. Combin- 
ing these considerations and our theoretical findings described above, we see 
that the T d 2 ^ 3 law is expected to be quite universal. 

2/3 

However, neither the T d law, nor even any other algebraic law were 
reported in the numerical-experimental studies in [HI E] • We think that the 
main reasons are that the algebraic fitting functions were not used, and, 
besides, solely the initial part of the distribution was built in [8], though the 
time range of the simulations allowed one to study the tails. Another point is 
that the tail of the disruption time distribution in the given problem should 
be considered separately from the initial part, because it corresponds to a 
different dynamical situation: here the regime of decay might be Poissonian 
(see analogues in [TjJJ 1201 EH]), or, in the very beginning, inverse square- root 

[MIESIITJ. 

Our theoretical inferences seem to be confirmed by the results of the very 
recent study [52], where the statistics of the decay process in the equal-mass 
three-body problem with randomized initial conditions were investigated in 
extensive numerical experiments. The lifetime distributions obtained in [52] 
have turned out to be heavy-tailed, i.e., the tails have turned out to be 
algebraic. The computed power-law index a for the integral distribution has 
been found to be within the narrow range, approximately from 0.4 to 0.7, 
depending on the virial coefficient (see [52]). The theoretically predicted 
value a = 2/3 is within this narrow range. 

The range of applicability of the derived "Tl — Td" relationship is similar 
to that of the derived distribution law. The proximity of our theoretical 
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quasilinear "Tl — TV' relationship to the numerical-experimental results in 
[S] (obtained for an equal-mass three-body system) has a natural heuristic 
explanation, consisting in that the stage before disruption is hierarchical, 
with the outer body exhibiting the usual final Levy flights, and this process 
can be described in the "hierarchical resonant scattering" setting for a three- 
body interaction (though, in general, triple encounters can play role). 

Finally, there is no wonder that a simple one-parameter two-dimensional 
map, such as the Kepler map, is able to describe the essential dynamics 
of disruption of a system with several degrees of freedom. The matter is 
that we consider a very special stage of the system evolution subject to se- 
rious limitations: the orbit of the escaping body is highly eccentric, and its 
pericenter distance is much greater than the size of the main binary. In 
many respects, reduction of the motion to the Kepler map in the problem 
considered is similar to reduction of the motion in the vicinity of the sep- 
aratrices of the "guiding resonance" to the separatrix map (see [12]) in a 
general Hamiltonian system. An important difference, however, is that the 
ordinary separatrix map has two parameters, and cannot be reduced, oppo- 
site to the case of the Kepler map, to a one-parameter form, because the 
phase increment in the ordinary separatrix map is logarithmic. 

7 Conclusions 

We have considered statistics of the disruption and Lyapunov times in an 
hierarchical restricted three-body problem. As we have seen, at the edge of 
escape the orbital periods of the escaping body exhibit Levy flights. Due to 
them, the distribution of the disruption times is heavy-tailed with the power- 
law index of the integral distribution equal to —2/3 , while the relation be- 
tween the Lyapunov and disruption times is quasilinear. The former finding 
is in accord with heuristic and numerical-experimental results in |JTJ W2\ H3]. 
while the latter one is in accord with recent numerical-experimental results 
in (HI [9] . Our theoretical results are valid for any system described by the 
Kepler map, i.e., they are valid at least in the external restricted three-body 
problem, both planar and non-planar, where the tertiary does not suffer close 
encounters with the central binary. The derived statistical laws appear to be 
valid as well in a more general "hierarchical resonant scattering" setting for 
a three-body interaction. 

The sharp difference between the two kinds of Hamiltonian intermittency, 
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in what concerns the slope indices of the asymptotic power-law distributions 
of the Poincare recurrences, allows one to explain the observed difference 
in the power-law indices of the distribution laws reported for the chaotic 
dynamics of the Solar system minor bodies. Dones et al. [39] reported on 
the algebraic tails of the integral distributions with the power-law index a 
equal to 0.8 ± 0.2, whereas Shevchenko and Scholl [THl E0] reported on the 
tails with the index equal to ~ 1.5. In [31] the escape times in the highly- 
eccentric chaotic cometary dynamics in the Solar system were studied (the 
perturbations due to the four giant planets were taken into account), whereas 
in [T9| [20] the subject was the low-eccentricity intervals between the eccen- 
tricity jumps in the chaotic asteroidal dynamics (in the restricted three-body 
problem Sun- Jupiter-asteroid). Judging by the values of the power-law in- 
dex, the former statistics correspond to the Hamiltonian intermittency of 
the first kind, whereas the latter one to that of the second kind. The statis- 
tics are LFl-dominated and LF2-dominated, respectively, and the predicted 
power-law indices a are respectively equal to 2/3 and « 3/2. The evident 
"inverse symmetry" of the indices is a property of the gravitational dynam- 
ics, described by the Kepler map; generally, when a separatrix map with an 
arbitrary 7 is in action, the predicted indices are I/7 and m 3/2, respectively. 
Thus the latter index is more "universal" than the former one. 

The change of the power-law index value (« 1 instead of ~ 2) in the 
"Tl - - Tj" relationship in comparison with the results on the dynamics 
of the Solar system minor bodies obtained in [TJEJEJIIJEJEJII] i s due 
to the fact that the singularity at crossing the separatrix is much stronger 
in the considered problem (which is described by the Kepler map instead 
of the ordinary separatrix map), and therefore the first kind Hamiltonian 
intermittency dominates over its second kind and thus defines the properties 
of the "Tl — T^" relationship and the tail of the disruption time distribution. 

We have shown that a critical non-zero value 7 cr it of the 7 parameter of 
general map fll]) exists that separates the maps with LF2-dominated dynam- 
ics from those with LFl-dominated dynamics. 

As a by-product of our study we have proposed a new method for estimat- 
ing the critical exponent a in Eq. (J7J). This method is based on computation 
of the 7 transition value separating maps (j3J) with LFl-dominated and LF2- 
dominated statistics of long recurrences. 

Finally, our inferences shed light on the mechanism of disruption of an 
hierarchical three-body system in the adopted setting of the problem: the 
typical way of disruption, as described by the Kepler map, is a kind of "Levy 
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unfolding" of the system in both time and space: at the edge of the system's 
disruption, the escaping body exhibits Levy flights in its orbital period and 
semi-major axis, and in the course of this random process the orbital period 
and semi- major axis become arbitrarily large until the separatrix separating 
the bound and unbound states of the motion is crossed and the body escapes. 
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Figure 1: A fragment of a chaotic trajectory of the Kepler map, demonstrat- 
ing Levy flights of the first kind (the prominent peaks of log 10 P or b on the 
left) and a Levy flight of the second kind (the oscillatory "plateau" of y and 
lo 

Sio -^orb cm the right), y is shown with a thin line, and log 10 P or b with a 
bold one; A = 5; the value of io is some big number chosen to exhibit this 
part of trajectory. 
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Figure 2: Integral distributions of the recurrence times measured in map 
time units (iterations) for the Kepler map; A = 5 (bold line) and A = 70 
(thin line). The straight dotted line shows the T~ 3 / 2 law. 
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Figure 3: The same as in Fig. |2j but the recurrence times are measured in 
real time units. The bold line corresponds to A = 5, and the thin one to 
A = 70. The straight dotted line shows the T" 2 / 3 law. 
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Figure 4: Integral distributions of the recurrence times measured in real time 
(bold lines) and map time (thin lines) units for the L-map (a) and for the 
Fermi map (b); A = 5; n it = 10 11 . The straight dotted lines show the T~ 3//2 
and T" 1 laws. 
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Figure 5: The dependence of the computed power-law index a on the 7 
parameter, for c = c reg (bold line) and c = (thin line); A = 10. 
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Figure 7: A statistical relation "log 10 Tl — log 10 ^r" , where Tl and T r are 
measured in map time units; A = 5. The straight dotted line shows the 
quadratic dependence. 
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Figure 8: The same as in Fig. [7J but the recurrence times are measured in 
real time units. The straight dotted line shows the linear dependence. 
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Figure 9: A spectrum of winding numbers, visualizing the resonant structure 
of the chaotic motion near chaos border; A = 5. 
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